Dynamics of the Bose-Hubbard model: transition from Mott insulator to superfluid 
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We study the dynamics of phase transitions in the one dimensional Bose-Hubbard model. To 
drive the system from Mott insulator to superfluid phase, we change the tunneling frequency at 
a finite rate. We investigate the build up of correlations during fast and slow transitions using 
variational wave functions, dynamical Bogoliubov theory, Kibble-Zurek mechanism, and numerical 
simulations. We show that time-dependent correlations satisfy characteristic scaling relations that 
can be measured in optical lattices filled with cold atoms. 



I. INTRODUCTION 

The spectacular experimental realization of the Bose- 
Hubbard model (BHM) using cold atoms in an opti- 
cal lattice |l| triggered an avalanche of both theoretical 
and experimental activity 0,13 • The excitement comes 
mostly from the fact that the derivation of the BHM in 
this system can be carried out rigorously @, Q, its pa- 
rameters can be experimentally manipulated in real time 
0, and lattice geometry can be engineered almost at 
will: it can be one, two, three dimensional, and can have 
different shapes, e.g., rectangular, triangular, etc. 

Physics of the Bose-Hubbard model is of both funda- 
mental and practical interest. Indeed, the BHM is one of 
the model systems on which our understanding of quan- 
tum phase transitions (QPTs) is based [HQ. The quan- 
tum phase transition happens in the BHM between the 
gapless superfluid (SF) phase and the gapped Mott in- 
sulator (MI) phase. Recently its signatures have been 
experimentally observed Q. In a homogeneous system 
at fixed density, the transition takes place only when the 
number of atoms is commensurate with the number of 
lattice sites. The practical interest in the BHM orig- 
inates from the possibility of realization of a quantum 
computer in a system of cold atoms placed in an optical 
lattice 0. 

In spite of experimental studies of the BHM and the 
large number of numerical and analytical contributions, 
understanding of the BHM physics is far from complete. 
In particular, a theory of the dynamics of the MI - 
SF quantum phase transition is still in its initial stages 
H, U 03 E3 ■ This is not surprising, as until very re- 
cently H, El E H Q E3 , QPTs have been studied 
as a purely equilibrium problem. The recent progress 
in dynamical studies has been obtained after applying 
the Kibble-Zurek mechanism (KZM) [TEII3, which was 
successful in accounting for non-equilibrium aspects of 
thermod yna mical phase transitions |l8j. to the quantum 
case [TllTlll5llT9ll20|. 

In this paper we investigate the dynamics of the one di- 
mensional (ID) BHM, focusing on two-point correlation 
functions. To describe their time dependence, we develop 
and use a variety of analytical approximations. We find 
that the two-point correlations satisfy simple character- 



istic scaling relations that should be experimentally mea- 
surable. Finally, we check the accuracy of our predictions 
with numerical simulations. 

Section [H] presents the model and defines the quan- 
tities of interest. In Section IIIII we discuss predictions 
coming from a toy two-site model. Section IIVI l|V|) ana- 
lyzes scaling relations of correlation functions induced by 
fast (slow) changes of the tunneling coupling. 



II. THE MODEL 

In terms of dimcnsionless variables used throughout 
this paper, the Hamiltonian reads 
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J + h.c.) + 
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1), (1) 



where we additionally assume a density of one particle 
per site. Such a model should be experimentally acces- 
sible in a ring-shaped optical lattice [2l|, where the ge- 
ometry of the problem imposes periodic boundary con- 
ditions on QJ. Another setup for investigations of the 
Bose-Hubbard model JTJ will be prov ided by the ongoing 
experiment in the Raizen group |22j, where a ID homo- 
geneous model with open boundary conditions will be 
realized. Below, we will assume periodic boundary con- 
ditions for the sake of convenience, but in a realistic ex- 
perimental situation with a few tens of lattice sites the 
dynamics should be unaffected by the boundary condi- 
tions. 

The Hamiltonian is driven from the MI to the SF 
regime by a linear ramp of the tunneling coupling 



J(t) = -, 



(2) 



where tq is the quench time-scale [l7[ l23| . The evolu- 
tion starts at t = from the ground state of QJ, i.e., 
1 1 , 1 , . . . ) , where the numbers denote boson on-site oc- 
cupations. The evolution stops at t — TQj max , where 
J max 3> 1< Therefore, the system ends up very far away 
from the critical point separating MI and SF phases: 
J w 0.29 24]. Experimentally, the change of the tun- 
neling coupling alone can be achieved by proper manipu- 
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lation of the lattice potential amplitude, followed by ad- 
justment of the atomic interaction strength via Feshbach 
resonances pHj . 

We are interested in the correlation functions: 

Ci{t) = ^{i>(t)\a\ +l a i + h.c.\rl>{t)), 

which are directly experimentally measurable because 
momentum distribution of atoms in a lattice is their 
Fourier transform ~ e~xp(ikl)Ci (k is the atomic mo- 
mentum) |2( . This observation shows that the correlation 
functions are good observables for our problem: by the 
end of time evolution J S> 1 so that interactions between 
atoms are asymptotically negligible. As a result, the cor- 
relation functions take well defined final values. 



III. DYNAMICS OF TWO SITE 
BOSE-HUBBARD MODEL 

In this section we consider a toy 2-site model, a prob- 
lem that can be completely solved analytically. The re- 
sults of this section will be useful later for studies of larger 
systems. Using symmetries of the Hamiltonian, one can 
show that the evolution starting from the uniform "Mott" 
state 11,1) leads to 



m))=a(t)\l,l)+b(t) 



where |a| 2 + \b\ 2 = 1 and 



|2,0) + |0,2) 
V2 



d 



dt 

A change of basis 




(a',b') 



(a-b,-a-b)/V2 



yields 



dt 
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19. 

4 ' 



(3) 



(4) 



(5) 
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This is exactly the Landau-Zener (LZ) model |2f|, whose 
relevance for dynamics of QPTs was recently shown in 
Refs. Q E [H HI- The quantity of interest is 
Ci(t) = 2\b'(t)\ 2 — 1, where b'{t) is provided by the exact 
solution of the Landau-Zener model in the case when the 
system starts its time evolution from the ground state at 
t = 0, i.e., from the anti-crossing center l20l|. This 
solution is a superposition of Weber functions (see Ap- 
pendix B of Ref . 20] ) , and it leads to 



Ci(oo) = -1 



sinh 



(tK" ,! 



r i 



(7) 



which has the following small tq (fast quench) expansion 



(8) 



For large tq (slow quench), we expand the gamma func- 
tions for large absolute values of the argument [2j| , 



T(z) = V2^z z ~h z ( 1 + — + — + Oiz 
y ' v \ 12z 288z 2 V 



and use that 



|r(zx)| 2 = 



r i - + i.v 



1 + ix 



IX 



x sinh(7ra;) ' 

7T 

cosh(7rx) ' 
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to obtain 



Ci(oo) 



+ 0{tq A ). 



(9) 



Eq. © is surprising since Ci(oo) = 1 — 2p ex , with p ex 
the excitation probability of the LZ system © at t = oo. 
Indeed, it implies that the excitation probability equals 



Pex(t : 



4 

72" 



when the LZ system (0 starts evolution from anticrossing 
center (t = 0) and evolves slowly till t — oo, while it is 
exponentially small (assuming tq ^> 1) 



p ex (t : -oo — > oo) = exp f- 



8 J 



for the LZ model evolving from t — — ooto< = oo. 

We have verified numerically that a slightly larger sys- 
tem (4 atoms in 4 lattice sites) exhibits the same scaling 
of Ci(+oo) in the fast and slow transition limit. Thus, 
these characteristics are not specific to a 2-site toy sys- 
tem only. In the following sections we will use different 
techniques to argue that the same scaling properties are 
shared by large lattice models. 

Before proceeding further, however, we mention that 
the power-law behavior of the excitation probability 
when the evolution starts at the anticrossing can be rel- 
evant for quantum adiabatic alg orithms psj . Indeed, by 
starting (or, by symmetry 20], ending) the algorithm 
near the anticrossing center the computation has a much 
higher failure probability. Thus, such situations have to 
be fiercely avoided when designing a path in Hamiltonian 
space between the initial and the solution Hamiltonians. 
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FIG. 1: Scaling properties of the first correlation function 
obtained numerically. Solid line: tq — 0.001, dots: tq — 0.03. 
Inset: solid line is a power law fit to data for 0.001 < tq < 0.1 
giving Ci(oo) = 0.501r£ 498 . All data is for M = 10 and 

Jmax — 600. 



In the following we will explain the observed behavior 
of C\ first by the time-dependent perturbation theory, 
and then by developing a bosonic Bogoliubov theory. 



A. Short time diabatic dynamics 

For short times we can approximate the wave function 
using time-dependent perturbation theory, 

m)) = o(t)|l, 1, • ■ •) + Kt) (|0,2, 1, . . .) + |2,0, 1, . . .) 
+ |1,0,2,1,...) + |1,2,0,1,...> +---)/V2M, (12) 

where M > 2 is assumed and |a| 2 + |6| 2 = 1. A time- 
dependent variational principle predicts in this case that 
dynamics of a(t),b(t) is governed by Eq. Q with tq 
replaced by tq/VM. Therefore, the familiar LZ problem 
shows up again, and it is useful to define new amplitudes 
a' and b' in the same way as in (JjJJ. Dynamics of a'(t) 
and b'(t) is governed by Eq. © with r — ► r/yM. 
To describe the build up of 



IV. FAST TRANSITIONS 

In this section we consider systems undergoing fast 
(tq <C 1) quenches. Let us start by summarizing some 
relevant numerical findings on C\. We studied numer- 
ically system sizes M = 3 • • • 10 (M is the number of 
lattice sites/atoms), and found that in all cases [2{| 



Ci(oo) 



OtT r 



(10) 



for tq's smaller than about 10 _1 . Depending on the sys- 
tem size, a € (0.37, 0.5) while j3 equals 1/2 within fitting 
errors: see the inset of Fig. ^for the M — 10 case. 

Moreover, as depicted in Fig. ^ the whole C\{t) func- 
tion after the rescaling 



Ci(t) 



(11) 



C 1 (t) = (2\b'(t)\ 2 -1)/VM 

for the wave- function l|12|) , we expand the exact solution 
of b'(t) |2fJ for small tq obtaining, in the lowest order, 



(13) 



Expression l|13|l is interesting: it implies that the way in 
which the first correlation builds up over time is indepen- 
dent of system size and takes some universal (indepen- 
dent of tq) form after simple rescalings. In Fig. [3 this 
prediction is compared to the numerical solution of the 
10-site Hubbard model. A perfect agreement is found for 
times smaller than about \^/tq. As will be explained 
in Sec. IIV Bl the number fluctuations start to develop 
significantly around t ~ -JjQi so it is not surprising that 
a simple wave function (|12|) fails to describe subsequent 
dynamics. 



takes an universal form for tq smaller than about 10 _1 . 

Two remarks are in order now. First, the two site 
prediction, Eq. |JSj , shares the same scaling with tq and 
a prefactor of the same order of magnitude (^/tt/A sw 
0.44) as the numerics for larger systems. Second, it is 
interesting to ask weather the scaling relation (|10J) can 
be experimentally verified. Taking 10 _1 as the largest 
tq for which i|10[l works very well, we get Ci(oo) w 0.16 
in M = 10 case (Fig. This is to be compared to the 
ground state predictions at (i) the critical point (Ci w 0.8 
[30|). and (ii) the asymptotic value deep in superfluid 
(Ci = 1). Thus, our results suggest that, despite the 
fast drive of the system through the transition point, the 
first correlation builds up macroscopically. Therefore it 
should be experimentally measurable. 



B. Bogoliubov theory 

Using the insight gained from the above studies, we de- 
velop Bogoliubov approach that includes a macroscopic 
number of excitations into the wave function and is able 
to describe longer than nearest neighbor correlations. 
Our approach is a variant of the theory developed by 
Altman and Auerbach [ll| for large density of particles. 

We truncate the Hilbert space to states with only 
{0, 1,2} particles per site. The initial state is the Mott 
state with exactly 1 particle per site. In a fast transition 
(tq < 1) we can get well into the superfluid regime of 
J 1 before any substantial number fluctuations have a 
chance to build up around the initial Mott state. Thus, 
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FIG. 2: Short time dynamics of the first correlation function. 
Numerics for M — 10 and J ma x = 600 is given by solid line 
(jQ = 0.001) and large dots (tq — 0.1). The dashed line 
presents Eq. 11131 . The pluses (+) stand for a numerical so- 
lution in the Bogoliubov model for tq — 0.001 and prediction 
lHH for I = (both data overlap). 



in a fast transition the truncation remains valid well in 
the superfluid regime. 

As already mentioned, the correlators C; are conserved 
by the hopping term in the Hamiltonian. The hopping 
term dominates when J 3> 1 and this is why in this 
regime the correlators are observed to be more or less 
constant, see Fig. ^ Our idea is to use a truncated 
theory to predict correlators Ci(t) up to an instant t so 
large that J it) ^> 1, but small enough to keep the num- 
ber fluctuations negligible. The predicted correlators do 
not change in the following evolution dominated by the 
hopping term, so that C;(f) « C;(oo). 

In the truncated Hilbert space we call 2 particles 
in a site a quasiparticle, and an empty site is called 
a quasihole. The Mott state with 1 particle in each 
site is now the "empty" vacuum state. We introduce 
the quasiparticle and quasihole creation operators as c\ 
and d\, respectively. Their action is best illustrated by 
mapping the boson occupation number onto two num- 
bers, (n c ,nd), where n c {rid) are quasiparticle (quasi- 
hole) occupation numbers. This way we have in each 
site: |2) = |(1, 0)>, |1) = |(0,0)), and |0) = | (0, 1)> . Since 
within the {0, 1,2} subspace we cannot have two quasi- 
particles/holes in the same site, the hard-core constraint 
has to be implemented: (c\) 2 = and (d\) 2 = 0. All this 
leads to Cj|..., (n c ,nd)i, —) = S lno \..., (n c - l,rid)i, ...), 
cf|— , (n c ,n d )i, ...) = 5 n c \—, (n c + •■■), and ana- 

logical relations for the action of di and d\ operators. Ad- 
ditionally, we have to remove the states with one quasi- 
particle and quasihole in the same site, |..., (1, ...), 
since these states also do not map onto the {0, 1, 2} sub- 
space. This is done by the projector P = \\i(^~ cfcidj di). 
Finally, we note that since we deal with hard-core bosons 



all the commutators of quasiparticle/hole operators at 
different lattice sites commute. 

In this new language the Hamiltonian Q in the 
{0,1,2} subspace equals exactly PH2P, where H2 is 
quadratic in quasihole/quasiparticle operators 



Ho = -J 



E [ 2 &i 



d\d, 



V2(diCj +d\c ] 



1 3 



E 



(14) 



where (i, j) denotes nearest neighbor pairs. We also men- 
tion here that the new operators satisfy periodic bound- 
ary conditions as the original bosonic operators do. 

The truncated Hamiltonian PH2P is exact in the 
{0,1,2} subspace, but it is not quadratic in c and d. 
In order to proceed we approximate H w H2 and lift 
the hard-core bosonic constraint in all subsequent cal- 
culations: from now on [ci,cj] = <5y and = Sij. 
This way we arrive at a bosonic theory with a quadratic 
Hamiltonian H2 leading to solvable linearized equations 
of motion. The quadratic theory remains self-consistent 
as long as average density of excitations [3lJ 



(did,) « 



1 



(15) 



When p cx > i it is likely to find quasiparticles and quasi- 
holes occupying the same lattice site and the constraint 
imposed by the projector P must be violated. 
We proceed by going to the momentum space 



1 



<M 



E 



c k e 



d r = 



1 



'M 



fee 



ikr 



To simplify time-dependent calculations we add to 
the Fourier transformed Hamiltonian two terms, 

-cj,Cfc), that both 



J2k J{c\c k -d\d k ) cos A; and J^k 
commute with the Hamiltonian itself and therefore do 
not change the evolutions considered in this paper. The 
resulting Hamiltonian becomes 

H 2 = — JE C0S ^ 3c£,Cfc + 3d\dk + 2v / 2(c J t<i_ /t , + h.c.) 



k 

It can be conveniently rewritten to the form 

H2 = -<*-*] 



(16) 



h — 3 J cos k — 2\/2Jcosfc 
— 1 

2 



2V2J cos/c 3Jcosfc 





Ck 




. dU . 



( 3 J cos k 



(17) 



Below we look at description of time evolutions leav- 
ing the discussion of static properties of our theory to 
Appendix [S] As in former sections we start time evolu- 
tion from J — 0, and J(t) is given by We work in 
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the Heisenberg picture where the system wave function 
(the ground state at J = 0: |1, 1, ...)) is time-independent 
while the operators evolve according to i-^c k — [c k ,H 2 ] 
[dk, H 2 \. It leads to 



and i-jjdk 



i — 

dt 







Ck 




. dl k _ 





1 — 3— cos k -2^2— cosfc 

2 TQ TQ 

2^2— cos A: 3— cos/c — i 

TQ TQ 2 



Ck 

dlk 



(18) 



which has the following general solution: 

dl k {t) = v k (t)c k (0)+u k (t)S_ k (0), 

with initial conditions u k (0) = u k (0) = 1, v k (0) = 
v k (0) = 0. After some algebra based on l)18|l one finds 
that v k = v k and u k = u* k - this simplification showed up 
thanks to the convenient addition of the two constants 
of motion to the Fourier transformed quadratic Hamil- 
tonian (see above). The time evolution of the modes is 
given by 



dt 



u k 









I - 3— cos k 

l TQ 

2y/2— cos k 

TQ 



-2-v/2— cosfc 

TQ 

3— cos/c — i 

r Q 2 



(19) 

Additionally, we see that the Bose commutation between 
the time-dependent operators requires that |wfc(t)| 2 — 
|i>fc(i)| 2 — 1, which is conserved by the time evolution 
(|19fl . All expectation values can be calculated after solv- 
ing (|19fl using the fact that the wave-function in the 
Heisenberg picture is \9) = |1, 1, ...) for all times, so that 
d r (0)|*) = and c r (0)|#) = 0. 

In the following we use perturbative solution of Eqs. 
(|19fl in powers of ^Jtq. The discussion is simplified by 
introducing a new time-like variable 



(20) 



whose form is motivated by the scaling property (jllfl . 
Equations l(T§|) become 



ds 



u k 

Vk 



cos k 



I V2 

-V2 -I 



Uk 

t'fe 



4Vi 



" 1 




u k 


-1 




Vk 



(21) 



Therefore, p ex <C 5 for s <C so that the 



for small s 

quadratic approximation breaks down at 



1 

71 



(22) 



or at t ~ \J t qI V 7 ^- In a linear quench © this break- 
down time corresponds to 

1 



J : 



» 1 



V2t q 



which is well in the superfluid regime for a fast transi- 
tion. Therefore, when tq <C 1, our linearized Bogoliubov 
approach does not break down until well in the superfluid 
regime. 

These calculations prove that Bogoliubov approach 
works reliably before J ^> 1 and the correlation func- 
tions are (see Appendix lAl for static predictions) [12 



Ci 



dk 
27 



(kl) h\v k \ 2 + V2{u k v* k + u* k v k )} . (23) 



Solving JUJ) we find that C 2t (t) = 0(tq), while 



C 2l+ i{t) _ Svr s 3 / 2 



OO 

E 

n=l 



„2l! 



(24) 



with coefficients 



OLln 



T(2 + 2n)r [1-21 + 2nr 1 T [3 + 21 + 2n\ 



n] r [i + r 



(x)„ = r[.x + n]/r[a;]. To obtain this series expansion 
we differentiate ^C/ in Eq. remove the resulting 

s-derivatives with the help of equations of motion l|21|l . 
keep only the leading terms ~ ^Jtq, and finally integrate 
such obtained 4-G\ over s to get C;(s). 

The first term in the I = version of (|24|l reproduces 
Eq. (jT3j) . As shown in Fig. [3 Eq. (|2"4"|) works perfectly 
until s ~ l/\/2, i.e., up to the expected breakdown of the 
Bogoliubov approach (|2*2*|l . 

Since we consider tq <C 1, we have J = 1/tq ^> 1, 
and thus the rest of evolution after J is dominated by 
the hopping term that does not change the correlation 
functions. Therefore, the correlators at the break-down 
time t are good estimates of the final correlation function: 



with Ufc(0) = 1 and v k (0) = 0. These equations can 
be solved iteratively in powers of the small parameter 

As a self-consistency check, we calculate the density 
of excitations, Eq. (|15fl . Assuming fast transition limit, 
tq <C 1, we solve Eqs. <|21[) to zero order in ^Jtq and find 
that 



1 



= — / dk\v k \ 



Cj(oo) « Ci(t). 



(25) 



Setting s = s = l/\f2 in 124(1 for definiteness we get 
with accuracy of 0(tq 2 ) 



Ci(oo) « 3.9 x lO- 1 ^, 
C 3 (oo) w -3.5 x 10- 3 ^, 
C 5 (oo) w 1.4 x 10~ 5 
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By solving Eas. (|21fl perturbatively we find with the ac- 
curacy 0(tq) that 

C 2 (oo) « 4.0 x 10- 2 t q , 
C 4 (oo) w -3.2 x 1CT 4 tq, 
C* 6 (oo) « 1.1 x 1(T 6 tq. 

The first correlation Ci(oo) fits well our numerical re- 
sults - compare to (|10|l . Reliable numerical verification 
of longer range correlations would require calculations 
done on systems larger than our M < 10. Indeed, in the 
small size numerics it is hard to filter out finite size ef- 
fects especially when the long range correlations, which 
are small in magnitude, are considered. Nevertheless, the 
Bogoliubov theory and our finite size numerics agree that 
correlations Ci decay fast with the distance I. 

In contrast to the simple 2-site toy model of Sec. 
IIIII the Bogoliubov approach is able to describe not 
only nearest-neighbor but also longer range correlations. 
However, it turns out that in fast transitions the corre- 
lation functions are dominated by the nearest-neighbor 
term C\ with other terms being relatively small, if not 
negligible. This explains why already the simple 2-site 
toy model gives such surprisingly accurate predictions 
for larger systems in the fast transition limit. 

V. SLOW TRANSITIONS 

In this section we focus on the limit of slow transi- 
tions, i.e., tq > 1. Numerical studies in this regime are 
extremely time consuming, therefore we concentrate only 
on analytical results. Our predictions are based on the 
Kibble-Zurek mechanism that was successful in describ- 
ing non-equilibrium thermodynamical phase transitions, 
and apparently works for quantum phase transitions as 

well El III EE EIH. 

According to KZM, excitations of the system after a 
slow transition have the characteristic length-scale |2^| 



(26) 



where z and v are critical exponents and the quench 
time tq is taken as {dJ/dt) -1 at the critical point (J 
is the parameter driving the transition). For the Bose- 
Hubbard model the dynamical exponent z = 1. The MI 
- SF transition (at fixed integer density of atoms) in a 
d-dimensional Bose-Hubbard model belongs to the uni- 
versality class of a (d + l)-dimensional XY spin model 
0. In one dimension this mapping implies that v — > oo 
(Kosterlitz-Thouless transition). As a result, 1 — G\, 
which is proportional to the hopping energy of long wave- 
length excitations, should scale for tq > 1 as 



i-ci(oo)~r 



(27) 



Q 



The exponent —2 means a rather steep dependence of 
the hopping energy on the quench time tq , which should 
make it easily discernible experimentally. 



Using i|26|) and (|27|l it is easy to provide predictions for 
two iv ~ 0.67 33]) and three (y = 1/2 (6J) dimensional 
Bose-Hubbard models. In the two dimensional case one 
has 



l-Ci(oo) 



1 



while in the three dimensional model 



l-d(oo) 



2/3 ' 



It would be very interesting to verify scaling relations 
shown in this section either experimentally or numeri- 
cally. 



VI. SUMMARY 

We described build-up of correlations in the BHM dur- 
ing transitions from Mott insulator to superfluid regime 
using: a variational wave function, the dynamical one 
dimensional Bogoliubov theory, the Kibble-Zurek mech- 
anism, and numerical simulations. The time-dependent 
correlations satisfy characteristic scaling relations that 
are directly experimentally measurable. 

This research was supported by US Department of En- 
ergy and NSA. J.D. was supported in part by Polish Gov- 
ernment scientific funds (2005-2008). 



APPENDIX A: GROUND STATE PROPERTIES 
OF BOSE-HUBBARD MODEL WELL IN THE 
MOTT REGIME 

Here we discuss the ground state properties of the 
Bose-Hubbard model predicted by the Bogolubov the- 
ory. The Hamiltonian l|16f) can be diagonalized by the 
Bogoliubov transformation 



c fc = u k B k + v* k Al k , d)_ k = v k B k 



u* k A J _ k , 



(Al) 



where u k (J) and v k (J) determine static properties of the 
Bogolubov vacuum. 

Here the Bogoliubov modes {u k ,v k ) are the eigenmodes 
of the stationary Bogoliubov-de Gennes equations 



u k 









-3 J cos k 



-2\/2Jcosfc 
l 

2 





u k 







(A2) 



2V2J cos k 3 J cos k 
with positive norm: \u k \ 2 — \v k \ 2 — 1 and eigenfrequency 



= - \J 4 J 2 cos 2 A; - 12 J cos fc + 1 . 



The normalization condition guarantees proper, i.e., 
bosonic commutation relations of A k and B k operators: 
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[i fe ,it] = [B k ,B}] = 6 kp , [A/i, B p ] = 0, etc. The diago- 
nalizcd Hamiltonian 

H 2 = «>* ( A l A k + B l B k) +J2 ( wfe + 3jc ° sfc - \ ) 

(A3) 

is a sum of Bogoliubov quasiparticle excitations. Its 
ground state is a Bogoliubov vacuum annihilated by all 
A k and B k . 

Now we calculate different quantities assuming that 
the system size M — > oo. As a self-consistency check we 
calculate the density of excitations 

Pox- 7^/ dk\v k \ 2 = 4J 2 + 0(J 4 ). (A4) 
The density remains <C 5 for J <C t^tj - compare to 

CE). 

The expression for correlation functions in the static 
calculations is obtained after using . Due to similarity 
in the notation, it is the same as (|23|l . except that now 
u k and v k depend on J rather than t. As a result we get 





= 4J + 0(J 3 ), 


c 2 


= 18J 2 + 0(J 4 ), 


c 3 


= 88J 3 + 0(J 5 ), 


c 4 


= 450J 4 + O(J 6 ), 


c 5 


= 2364J 5 + 0(J 7 ), 


Gs 


= 12642J 6 + 0(J 8 ), 


c 7 


= 68464J 7 + 0(J 9 ). 



These results in the lowest nontrivial order listed above 
agree perfectly with perturbative ones. Indeed, C\ ■ ■ ■ C3 
can be found in Eq. (5) of [3Cj , while C4 • • • C7 match 
unpublished results of one of us (B.D.). Discussion of 
this intriguing finding is beyond the scope of this study 
and is left out for further detailed investigations. 



To close discussion on static properties of our theory 
we notice that also the ground state enenergy per site (£), 
predicted by the Bogolubov theory, agrees in the lowest 
order with exact perturbative calculation. Indeed, we get 
from (|A"3l that 



£ = ^-l dk (u3 k + ZJcoak- M = -4J 2 + 0(J 4 ), 



which has to be compared to Eq. (3) of pjof. 
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